home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / vslesing.cc < prev    next >
C/C++ Source or Header  |  1995-12-20  |  2KB  |  75 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *    Verify application of SVD to solving sets of simultaneous
  6.  *                linear equations
  7.  *
  8.  * $Id$
  9.  *
  10.  ************************************************************************
  11.  */
  12.  
  13. #include "svd.h"
  14. #include <iostream.h>
  15.  
  16.                     // Solve Ax=b and test the result
  17. static void test_sle(const Matrix& A, const Vector& b, const Vector& x_true)
  18. {
  19.   SVD svd(A);
  20.   cout << "\ncondition number " << svd.q_cond_number() << endl;
  21.   Vector x = SVD_inv_mult(svd,b);
  22.   Matrix solutions(2,x_true.q_upb());
  23.   MatrixRow(solutions,1) = x_true;
  24.   MatrixRow(solutions,2) = x;
  25.   solutions.print("true vs. computed solution");
  26.   Vector b_comp = x;
  27.   b_comp *= A;
  28.   cout << "\tchecking to see that Ax is indeed b\n";
  29.   verify_matrix_identity(b_comp,b);
  30. }
  31.  
  32. static void test1(const int neqs, const int nvars)
  33. {
  34.   cout << "\n\nChecking solution of a set of linear equations\n"
  35.           "with a Hilb+E of order " << neqs << "x" << nvars << endl;
  36.   Vector x(nvars);
  37.   Matrix A(neqs,nvars);
  38.   A.hilbert_matrix(); (MatrixDiag)A += 1;
  39.   struct fill : public ElementAction
  40.   {
  41.     void operation(REAL& element) { element = i; }
  42.   };
  43.   x.apply(fill());
  44.   Vector b = x;
  45.   b *= A;
  46.   test_sle(A,b,x);
  47. }
  48.  
  49. static void test2(const int dim)
  50. {
  51.   cout << "\n\nChecking solution of a set of linear equations\n"
  52.           "with a Hilbert matrix of order " << dim << endl;
  53.   Vector x(dim);
  54.   Matrix A(dim,dim);
  55.   A.hilbert_matrix();
  56.   struct fill : public ElementAction
  57.   {
  58.     void operation(REAL& element) { element = i; }
  59.   };
  60.   x.apply(fill());
  61.   Vector b = x;
  62.   b *= A;
  63.   test_sle(A,b,x);
  64. }
  65.  
  66.  
  67. main(void)
  68. {
  69.  cout << "\n\nVerify application of SVD to solving sets of simultaneous "
  70.          "equations" << endl;
  71.  test1(10,10);
  72.  test1(20,10);
  73.  test2(10);
  74. }
  75.